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Abstract. SU(N) quantum spin systems may be realized in a variety of physical 
systems including ultracold atoms in optical lattices. The study of such models also 
leads to insights into possible novel quantum phases and phase transitions of SU (2) 
spin models. Here wc use Gutzwiller projected fermionic variational wavefunctions to 
explore the phase diagram and correlation functions of SU (TV) quantum spin models in 
the self-conjugate representation. In one dimension, the ground state of the SU (4) spin 
chain with Heisenberg bilinear and biquadratic interactions is studied by examining 
instabilities of the Gutzwiller projected free fermion ground state to various broken 
symmetries. The variational phase diagram so obtained agrees well with exact results. 
The spin-spin and dimer-dimer correlation functions of the Gutzwiller projected free 
fermion state with N flavors of fermions are in good agreement with exact and 1 /N 
calculations for the critical points of SU(N) spin chains. In two dimensions, the phase 
diagram of the antiferromagnetic Heisenberg model on the square lattice is obtained 
by finding instabilities of the Gutzwiller projected 7r-flux state. In the absence of 
biquadratic interactions the model exhibits long range Neel order for N = 2 and 
4, and spin Peierls (columnar dimer) order for TV > 4. Upon including biquadratic 
interactions in the 577(4) model (with sign appropriate to a fermionic Hubbard model), 
the Neel order diminishes and eventually disappears, giving way to an extended valence 
bond crystal. In the case of the SU (6) model, the dimcrized ground state melts at 
sufficiently large biquadratic interaction yielding to a projected 7r-flux spin liquid phase 
which in turn undergoes a transition into an extended valence bond crystal at even 
larger biquadratic interaction. The spin correlations of the projected 7r-flux state at 
N = 4 are in good agreement with l/N calculations. We find that the state shows 
strongly enhanced dimer correlations, in qualitative agreement with recent theoretical 
predictions. We also compare our results with a recent quantum Monte Carlo study 
of the SU (4) Heisenberg model. 
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1. Introduction 

Quantum spin models provide a setting in which one can explore interesting strong 
correlation physics that arises from quantum fluctuations. Such fluctuations can be 
large enough, in certain cases, to melt any form of classical order leading to various 
exotic spin liquid states[U [2]. A class of spin liquids of particular interest are those 
that support gapless excitations. Apart from experimental questions regarding their 
possible existence [31 HI El El [7] , it is interesting to inquire under what circumstances 
such gapless spin liquids can arise in two or more dimensions. The Bethe ansatz solution 
of the spin-1/2 Heisenberg antiferromagnet chain is a well-known exact example of such 
a gapless spin-liquid in one spatial dimension, but much less is known about higher 
dimensions. This question was raised in the early days of high-T c theory, as Anderson's 
original proposal for a resonating valence bond (RVB) spin liquid had a Fermi surface 
of gapless spinon excitations [H [9]. 

One route to accessing the physics of some of these spin disordered states is by 
generalizing the usual SU(2) spin models to SU(N) models. As N increases from 2 
to larger values, quantum spin fluctuations are enhanced, weakening any spin order. 
Such models can accommodate several types of spontaneously broken symmetries: Spin 
order, spin dimerization, and charge-conjugation symmetry breaking[10]. We note that 
models of SU(N) quantum spins are not purely theoretical exercises: It may be possible 
to realize them with ultracold atoms in optical lattices[TTl [121 m quantum dot 
arrays[ll], or as special points in models with spin and orbital degrees of freedom [15], 
although the effects of SU (N) symmetry breaking terms need to be carefully examined 
in each case. 

Of the many possible representations of SU(N), we focus on a particular self- 
conjugate representation with N/2 fermions on each site, each with a different flavor 
due to the Pauli exclusion principlepSl [17]. In this representation, the generators of 
the SU (N) algebra may be expressed in terms of the fermion creation and annihilation 
operators as: Sg(i) = fifi/3 — \o~p- The constraint S"(i) = thus holds in the subspace 
with exactly N/2 fermions on each site, and consequently there are the correct number 
(N 2 — 1) of special unitary generators. The representation is called "self-conjugate" 
because upon making a particle-hole transformation the same representation, namely 
one with N/2 fermions, is obtained. All representations of SU (2), regardless of the total 
spin, are automatically self-conjugate, but only certain representations of SU(N > 2) 
are self-conjugate. An advantage of the self-conjugate representation is that it is easy 
to construct SU(N) invariant Hamiltonians that retain all of the symmetries of the 
underlying lattice, and hence mimic SU(2) models in this regard. For more details, 
including the Young tableau classification, see reference [10] . 

In the N — > oo limit saddle point solutions of the fermionic path integral are 
exact and the operator constraint S®(i) = constraint can be replaced by the much 
simpler mean-field constraint (S^(i)) = 0. These large- N SU(N) antiferromagnets 
cannot break global SU(N) spin symmetry; some possess ground states that do not 
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break any lattice symmetry and thus furnish mean field caricatures of a class of spin 
liquids. Going from N = oo down to the physical limit of N = 2 requires the inclusion 
of fluctuations about the mean field state, and the problem can be recast in the form 
of a strongly coupled gauge theory interacting with fermionic matter fields. Progress 
toward calculating the properties of such a field theory relies on a 1/N expansion, 
and some results have been obtained in this manner for one and two dimensional 
models [HI UHl 1201 1^11 [331 |SHl ISOl tHXl E2I though not without 

controversy. There is also good reason to be concerned about the reliability of such 
an expansion in the physically important SU (2) N is of course no longer large. 

DMRG calculations for SU(N) quantum antiferromagnets [Tl] and Hubbard models [33] 
provide some confirmation of analytical understanding of one- dimensional (ID) chains. 
In two dimensions (2D) there is a very interesting quantum Monte Carlo (QMC) study 
of SU(N) quantum antiferromagnets by Assaad[3l] that we discuss further in sectional 

A different approach to the study of SU(N) spin models also begins with mean 
field states that satisfy the average constraint (S"(i)) = 0, but then the operator 
constraint S"(i) = is implemented via an on-site Gutzwiller projection. The Gutzwiller 
projection operator forces the number of fermions to be precisely N/2 at each site. For 
N = 2, and denoting as usual the two spin states by f, j, the projection operator 
takes the well-known form Pq = ~ n i] n ii)- The resulting many-body wavefunction 
thus lives in the correct Hilbert space for SU(N) antiferromagnets and serves as a 
variational approximation to the spin ground state. For N = oo, a probabilistic 
central limit argument shows that the Gutzwiller projection is unimportant and gives 
the same result as the mean field state. For any finite N, projection is crucial 
and nontrivial; the advantage of the approach is that the projection constraint can 
be handled numerically exactly with the variational Monte Carlo algorithm. Thus 
the generalization to SU(N) provides a rationale for understanding the Gutzwiller 
procedure: Mean-field wavefunctions obtained in the large- iV limit are then modified 
by projection to account for the occupancy constraint at finite-iV. The approach 
however suffers from the criticism that it is biased and restricted by the choice of the 
variational mean field state (or equivalently the preprojected wavefunction) and that 
local Gutzwiller projection may not account for all of the important correlations. 

The Gutzwiller variational approach has been applied to a variety of SU{2) 
quantum antiferromagnets; see for instance Refs. [351 EE [37J [381 ESI HQl SU H21 03]. 
Many of these studies support the possibility of 2D spin liquids with gapless spin 
excitations. One advantage of extending the Gutzwiller variational approximation to 
SU(N) quantum antiferromagnets is that the approximation becomes more accurate 
as N increases. Indeed, since the gauge theory approach and the variational approach 
reduce to the same (exact) mean field theory at N = oo, but suffer different criticisms 
at finite N, it is interesting to compare results from both approaches as N is decreased 
systematically from N = oo down to N = 2. The existence of some exact results for 
SU(N) quantum antiferromagnets in one and higher spatial dimensions [10 j provides 
valuable checks not available in the SU (2) case. So motivated, we numerically explore 
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in this paper Gutzwiller wavefunctions for SU(N) spin models in the self-conjugate 
representation with Heisenberg bilinear and biquadratic interactions. In section [2] 
we compare the variational approach in ID with exact results and analytical 1/N 
calculations. Section [3] focuses on the 2D square lattice. Comparison is made with 
1/N calculations and with Assaad's quantum Monte Carlo results. We summarize and 
discuss the implications of our results in section HI For the reader interested in the main 
results, phase diagram figures in the different sections provide a quick overview of our 
conclusions. 

2. One Dimension 

We begin with a study of SU (N) spin models in ID. The existence of a reliable phase 
diagram for the SU(A) spin chain[1.0] provides a valuable check on the quality of the 
variational wavefunctions. We first define the model and then compare the phase 
diagram of the SU (4) chain as obtained with the variational wavefunctions to the known 
result. Then we examine various correlation functions calculated from the Gutzwiller 
projected SU(N) Fermi gas wavefunction and compare the exponents so obtained to 
exact analytical results for the critical point in the SU(N) spin chain. 

2.1. Models 

For N = 2, the usual Heisenberg model is the only nearest-neighbor Hamiltonian that 
is both SU(2) symmetric and translationally invariant. In ID we may write it as: 

H 2 = J2S^)S^ + 1), (1) 

i 

where the trace Tr[S(i)S(i+l)} is simply a rewriting of the vector form H 2 — 2 J2i Si-S i+ i 
with standard spin operators Si = \fi<j"f v i in terms of the matrix form of the 
spin operators. It is well-known that the model has no long range order, but rather 
exhibits algebraically decaying antiferromagnetic spin correlations up to a mutiplicative 
logarithmic factor. 

For N = 4, the most general translationally-invariant nearest-neighbor Hamiltonian 
can have an additional biquadratic spin-spin interaction term: 

H A = cos^E Sf(i)SZ(i + 1) + ^ Em) sP a(i + I)] 2 , (2) 

i i 

where —tt < 9 < n parametrizes the relative strength of the Heisenberg and 
biquadratic terms. Whereas the usual bilinear Heisenberg interaction exchanges two 
fermions on adjacent sites, the biquadratic term exchanges two pairs of fermions. The 
antiferromagnetic region of the phase diagram of H4 is generically gapped [TTl [TOj |4~4"] . 
The Lieb-Schultz-Mattis theorem then says that these gapped phases must break 
translational symmetry in one way or another [45] . Tuning 9 leads to a variety of 
phases and phase transitions. For instance, positive 9 frustrates dimerization, and when 
9 becomes sufficiently large, the dimerized phase is eliminated. 
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2.2. Phase Diagram of the SU(A) Model 

The phase diagram of the ST/(4) model, equation [2], was obtained in reference [10] and 
is shown in the inset to figure [U It displays four phases: A fully polarized ferromagnet 
(TAi), a dimerized phase (£>), a phase with broken charge-conjugation symmetry (C), 
and, finally, a phase with broken translational symmetry and a 6-site unit cell (6- fold). 
The T> to TM. transition and the 6-fold to TM. transition are both first order, while 
the T> to C and the C to 6-fold state transitions are both continuous. We now describe 
these phases more precisely and present the results of calculations based upon the 
corresponding variational wavef unctions. 

J-~M., Ferromagnet: The ground state of a fully polarized ferromagnet breaks the global 
SU(4) spin symmetry but no lattice symmetries. A simple and exact ground state 
wavefunction may then be constructed by placing any two of the four fermion flavors 
on the lattice, each site having the same two flavors. All other ground states can be 
obtained by global SU(4:) rotations of the state. Because the Pauli exclusion principle 
prevents any fermion hopping in the J-~A4 state, it is straightforward to show that 
(S%(i)S%(i + 1)) = ([S${i)SP(i + l)] 2 ) = 1, so that the exact ground state energy per 
site is 

e^f = cos0 + ~sin0. (3) 

J), Dimerized: The dimerized phase is a spin gapped phase with a two fold ground state 
and a nonzero order parameter 5-p = (S^(i)S^(i + 1)) — (Sg(i — l)S^(i)} that takes on 
equal positive/negative values in the two ground states and is zero in a state with no 
broken symmetries. This phase does not break the 577(4) spin symmetry, but does break 
lattice translations (the ground state is invariant only under translation by two lattice 
spacings), and inversion symmetry about a lattice site. In order to obtain a wavefunction 
for this phase, we consider a mean field fermion Hamiltonian with alternating hopping 
strengths (1 + 5) and (1 — 5) on successive bonds: 

H°v = ~ E (! + *(-!)*) [ftfi + i,a + h.c] . (4) 

i,a=l..A 

This starting Hamiltonian has the following favorable features: it is 5C/(4) symmetric, 
has a gap to fermion excitations (and thus a gap to spin excitations in mean field theory), 
and breaks the same lattice symmetries as the dimerized phase. In addition, the ground 
state of the Hamiltonian satisfies (S%(i))o = since it is particle-hole symmetric. The 
ground state of this mean field model is simply a product of four Slater determinants, 
one for each flavor of fermion, with the lowest half of the single particle states filled. 
Gutzwiller projecting the mean-field ground state leads to a variational ansatz for the 
dimerized phase of the spin model, with 5 being the variational parameter that we 
optimize by minimizing (H4) to find the best variational ground state. We find that the 
dimerization strength, defined in the mean field problem via the variational parameter 5, 
is nonzero at 9 = indicating that the S77(4) Heisenberg model has a dimerized ground 
state. With increasing 9, the optimal 5 decreases and vanishes around 9 ~ 0.41(3), in 
reasonably good agreement with the exact result 9 C = tan -1 (1/2) « 0.4636. We also 
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find that 5 increases in magnitude for negative 9 and approaches unity at large negative 
9 indicating that dimerization is nearly complete. 

C, Charge-conjugation symmetry broken: The model if 4 is invariant under the global 
particle-hole transformation f ia <-> f} a and thus possesses charge-conjugation (C) 
invariance. In terms of the spin operators, the transformation takes the form Sg(i) — > 
—S^(i). The C phase corresponds to a phase in which this symmetry is spontaneously 
broken. This phase does not break the SU(4) spin symmetry but does break lattice 
symmetries as the ground states are invariant only under translation by two sites, and 
they break inversion symmetry about the bond centers. It is characterized by an order 
parameter made up of a triple product of spins: (Tr(SiSi+iSi+2)} ■ An extreme caricature 
of the state (analogous to the product state of nearest neighbor dimers) is an extended 
valence bond solid of site-centered SU(4) singlets formed from two flavors of fermions at 
a central site combined with a fermion from each of the two flanking sites. This product 
state is an exact ground state of model if 4 at the special point 9 = 9*= tan -1 (2/3). 

The mean field Hamiltonian we use to obtain the preprojected wavefunction for the 
C-breaking state is 

H° c = - E [ftfi + i, a + h.c.]-t* J2 (-lY[f} a fi + 2, a + h.c] . (5) 

i,a=1...4 i,o=1...4 

Nonzero t* breaks the global particle-hole symmetry since it connects sites belonging to 
the same sublattice. The alternating sign of t* on the odd and even sublattices breaks 
inversion symmetry about bond centers of the lattice, and generates a gap in the mean 
field fermion spectrum (and thus a gap to spin excitations in the mean field theory). 
Finally, the Hamiltonian is invariant under a global particle hole transformation followed 
by translation by one lattice spacing, and hence satisfies (S%(i) + S%(i + l))o = 0. We 
can further modify the wavefunction to include a staggered chemical potential in order 
to obtain (S"(i))o = at each site. However since the mean field state is projected 
into the correct Hilbert space that satisfies S%(i) = exactly, we choose to work with 
the simpler mean field Hamiltonian without this staggered chemical potential (we have 
checked that including the staggered chemical potential does not affect the results in 
any significant quantitative manner). 

A check on the quality of this variational wavefunction for the C phase is provided 
by a comparison between the variational and exact ground state energy per site at the 
point 9* = tan~ 1 (2/3) where a C product is the exact ground state while the variational 
ground state exhibits nonzero if. The variational ground state energy at this point, 
E v3ur (9*) = -0.6934(10), is close to the exact result E CXSlCt (9*) = -^4- « -0.69337. . .. 

Turning to general 9, we find that the variational parameter t* is zero for 9 < 
0.42(2), and increases monotonically for 9 > 0.42(2). The 9 at which t* first becomes 
nonzero thus seems to coincide (within numerical error) with the point where the 
dimerization parameter 5 vanishes (9 = 0.41(3)). 

6-fold degenerate state: The phase diagram of the SU (4) spin chain has two special 
points with enlarged SU(6) symmetry, and a gapped 6-fold degenerate phase is 
associated with one of these pointspjO]. We may understand the origin of the 6- 
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Figure 1. Variational calculation of the ground state energies of dimcrized T> and 
charge-conjugation broken (C) phases of the model Hamiltonian H± in equation [2] 
Both broken symmetries apparently vanish around 9 = 0.42(3), hinting at a possible 
continuous phase transition with a gapless critical point that is reasonably well- 
described by the Gutzwillcr projected 4- flavor Fermi gas wavefunction (see text). 
Inset: The full phase diagram, including the ferromagnetic and 6-fold symmetry broken 
phases, of the SU(A) spin chain with biquadratic interactions. The exact location of 
the T> — C transition is at 9 = tan _1 (l/2) w 0.4636. The point at which the C extended 
valence bond product state is an exact ground state is 9 = tan -1 (2/3). 

fold degeneracy as follows: View the 6 possible states of the self-conjugate SU(4) 
representation (two distinct flavors chosen out of four possible ones) on each site as 
the six states of the fundamental representation of SU(6). Six such states, taken from 
six adjacent sites, can be combined into a SU(6) singlet. The resulting spin-gapped 
ground state breaks translational symmetry with a 6-site periodicity. Unfortunately 
we have not yet found a way to express the state in terms of the Slater determinants 
of single-particle wavefunctions, a proper description likely needs some form of pairing 
to capture the above physics. We therefore do not focus on this phase at present, 
postponing it to future study. 

T> — TM. transition We know from the phase diagram of reference [10] that there is 
a first order T> to TM. transition. Since the dimerization order parameter appears to 
increase monotonically at larger negative values of 8, we may estimate the approximate 
location of the T> — TM. transition by comparing the energy of the T M. state with 
a fully dimerized state (5 = 1) that is just a product of nearest neighbor singlets on 
alternate bonds. The energy of this variational state can be found knowing that for 
neighboring uncorrelated sites (not on the same singlet bond) (Sp(i)S^(i + 1)) = 
and {[S%(i)SP (i + l))] 2 ) = 5/3, while for neighboring spins that form the singlet bond 
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(Sp(i)S%(i + 1)) = -5 and ([S$(i)SP(i + l)] 2 ) = 25. This leads to 

4 =1 = --cos 9 + — sin9. (6) 

Comparing this energy with e^f, we find the variational estimate of the angle at which 
the dimer state becomes unstable to ferromagnetism to be 0™fm = tan~ 1 (84/74) « 
— 0.737T, very close to the exact result 0^ C ^ M = — 3n/A. The error in the variational 
estimate of the transition angle is consistent with the fact that while the energy of the 
ferromagnet is obtained exactly, the energy of the dimer product state with e^ 1 is only 
an upper bound to the true ground state energy of the dimerized state at the transition. 
C — T> transition The numerical coincidence of the values of 9 at which S and t* vanish 
indicates that the transition between the T> and C phases could be continuous even 
within the variational approach as in the rigorous phase diagram. We have not studied 
wavefunctions with coexisting t* and S broken symmetry parameters and cannot rule out 
the possibility that such a variational ansatz may have lower energy in the region close 
to the transition. We also cannot rule out the possibility that the variational approach 
leaves a very small window where t* = 5 = giving rise to a gapless phase. Assuming 
however that neither of these possibilities is realized, and that the transition between the 
two phases is continuous even within the variational approach, the projected half-filled 
Fermi gas state (with 5 = t* = in the mean field Hamiltonian) is a good candidate for 
the critical point describing the C — T> transition. We turn next to a study of correlation 
functions of this wavefunction. 



2.3. Correlation Functions for the Projected Fermi Gas at Various N 

As the phase diagrams of the SU{2) and SU(A) spin chains are known [TO], these provide 
good test cases for the method. At the critical point in the SU(N) chain, the exponents 
of the spin-spin and dimer-dimer correlation functions follow directly from equations 
(4.10) and (4.12) of reference [TO] , once the scaling dimension of the level k = 1 Wess- 
Zumino-Witten (WZW) field g is known. (The contribution of the currents Jl and Jr 
in equation (4.10) to the spin-spin correlation function is subleading, as the currents 
have dimension 1, greater than that of the g- field.) Reference [46] gives the dimension 
of a tower of WZW operators labeled by the integer a: 

dim(g) = h + h 

h =~ h = 5 ^,« = l,2,...,Ar, (7) 

with a = 1 corresponding to the operator with smallest nontrivial dimension. The 
conformal charge c = iV — 1 is correct as it equals the number of fermions in the 
corresponding SU(N) Hubbard model minus 1 due to the freezing out of charge 
fluctuations, leaving only N — l spin excitations. Now the exponent of the staggered part 
of the spin-spin correlation function, as well as the dimer-dimer correlation function, is 
twice the dimension of g, yielding 2 — 2/N. This exponent reduces to the free fermion 
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In (L) 

Figure 2. Logarithm of the staggered and uniform components of the spin-spin 
correlation function at distance L/2 plotted versus ln(L), for the projected free Fermi 
gas with N- flavors (L denotes the number of sites of the spin chain). The uniform spin 
correlations appear to decay as 1/r 2 for all N. The staggered spin correlations are 
enhanced at smaller N, and decay with the indicated exponents (see text for details). 



value of 2 in the N — ► oo limit, and to 1 in the usual SU(2) Heisenberg chain, as it 
should. 

How do correlation functions behave in the projected Fermi gas wavefunction for 
general N? For N = 2 the Gutzwiller projected Fermi gas wavefunction at half-filling is 
the exact ground state of the Heisenberg model with 1/r 2 interactions [47J HH]. It also 
correctly describes correlation functions of the ground state of the Jl — J2 Heisenberg 
model at the critical point between the gapless spin fluid phase and the dimerized phase 
[33]. The spin-spin correlations of this wavefunction have been computed exactly [50] [51], 
but we are not aware of an exact result for its dimer-dimer correlations. We present 
numerical results for both correlations below. For N = 4, given the possibility that the 
projected free Fermi gas wavefunction at half-filling could be a candidate for the C — T> 
transition, we study spin-spin and dimer-dimer correlation functions of the wavefunction 
and compare to exact results from the field theory for this transition. We also examine 
the correlation functions in the projected Fermi gas wavefunction for N > 4 as such 
wavefunctions may describe multicritical points in the phase diagram of generalized 
SU (N > 4) spin models. 

Spin-spin correlations: At N = oo, the long distance behavior of the spin correlation 
function C ss (x) = (Sg(0)S^(x)) is given by mean field theory, C ss (x) ~ (— 1/x 2 + 
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Css(x) ~ -— a . (8) 



(— l) x /x 2 ). In the opposite limit, N = 2, the spin correlations of the projected Fermi 
gas wavefunction have been calculated exactly by Gebhard and Vollhardt [501 EI] to be 
C ss {x) = (— 1)3 ^(ra) ; where S"i(x) is the sine integral function S'i(x) = Jq dy sin(xy) / y . 

The long distance decay of the correlator is C ss (x) ~ -jp. Thus the staggered 

spin correlations decay more slowly for N = 2 than they do in mean field theory, while 
the uniform component continues to decay as 1/x 2 . 

Incorporating gauge fluctuations perturbatively[20| |2T] modifies the long distance 
spin correlations to: 

A{-lf B 

The 0(1/N) result for the exponents are j3 s = 2 and a s = 2 — 2/N; rather surprisingly, 
the latter exponent agrees with the exact result. Thus, gauge fluctuations enhance the 
staggered spin correlations over the mean field result, while the uniform component 
decays at the same rate ~ 1/x 2 as the mean field result. 

In order to extract the exponent ot s (N) from the numerical calculations on the 
Gutzwiller projected wavefunction for general N, we assume that the correlations are 
of the form 

c-M = - § ■ 0) 

In order to carry out finite size scaling, we consider the function 

C(L/2, = (-l)^C„ (i /2, = ^-^y^. (10) 

From this equation we can obtain the staggered and uniform components of the spin 
correlations via 



Cstag(£) — 7^ 

C miii (L) = - 



C(L/2) + \ (C(L/2 + 1) + C(L/2 - 1)) 
C(L/2) - i (C(L/2 + 1) + C(L/2 - 1)) 



(11) 
(12) 



To leading orders in 1/L, these functions take the form 



r &t ( n ~ i v ; i (id 

stagl ; ~ (L/2) a * (L/2) 2 +^ (L/2) 4 1 ; 

unifl J ~ (L/2) 2 L 2 (L/2) a ° (L/2) 4 1 J 

We also obtain the staggered and uniform components of the spin correlation function 
data for various L using equation [12j We first fit the staggered correlations using the 
parameters A s ,u s ,a s . Next we use this value of a s and the parameters B s ,v s to fit to 
the uniform component of the spin correlations. The fits are shown in figure [2] for both 
the uniform and staggered components for cases iV = 2, 4, 6, and 8. This leads to the 
following estimates for various N: a s (2) = 0.996(4), a s (4) = 1.37(2), a s (6) = 1.63(2), 
and a s (8) = 1.74(2). These values are in remarkably good agreement with the exact 
value of the spin-spin correlation function exponent at the critical points of these chains, 
which for N = 4 lies at the continuous C — T> transition. We also find that the fitted 
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In (L) 

Figure 3. Logarithm of the alternating component of the dimcr-dimer correlation 
function at distance L/2 plotted versus m(L). for the projected free Fermi gas with 
N- flavors (L denotes the number of sites of the spin chain). The decay exponents 
are indicated (see text). The decay of the uniform correlation is much faster and the 
corresponding exponents have not been obtained reliably. 

value of the amplitude ratio A s /B s tends to unity as N — > oo, consistent with mean 
field theory. At present, we are unable to determine if the small differences between the 
exact results and the wavefunction calculations of a s (4) and a (6) are real or an artifact 
of working with chains of less than 100 sites. 

Dimer-dimer correlations: We have similarly evaluated the dimer-dimer correlations in 
the projected iV-flavor Fermi gas in ID, namely Cdd{x) = (Sp(Q)SP(l)Sf}(x)S"(x + 1)). 
The finite size scaling of this correlation function is analyzed in a manner similar to 
that of the spin-spin correlation function, except for one significant difference. We fit 
to Cdd(x) = Ad{— l) x /x ad , dropping the uniform component that decays much more 
rapidly, so that we cannot extract its behavior reliably compared to the staggered 
component. The finite size scaling plots of the correlation function are shown in figure 
[3j along with estimates of ad for various N. We find a^(2) = 1.02(5), a^(4) = 1.0(1), 
a<i(6) = 1-4(2), and a^(8) = 1.5(2). The projected Fermi gas wavefunction thus has 
strongly enhanced alternating dimer correlations in addition to enhanced staggered spin 
correlations. For N = 2, it is in agreement with the exact result ad = 2 — 2/N. For 
N > 2, the variational wavefunction exponents ad(N) < a s (N) (most significantly for 
N = 4) while exact results suggest a s (N) = atd(N) = 2 — 2/N. Thus, the projected 
Fermi gas wavefunction does not quite capture this aspect of the C — V critical point at 
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N = 4, although it does capture the existence of strongly enhanced staggered spin and 
alternating dimer correlations, decaying in power-law fashion with anomalous exponents. 
Having shown that the variational wavefunctions provide a reasonably good description 
of SU(N) spin models in ID, we next turn to 2D examples. 

3. Two Dimensional Square Lattice 

Buoyed by the successful description of the SU(N) spin chains with variational 
wavefunctions, we now apply the same methodology to the study of 2D SU(N) spin 
models. Much less is known reliably about the phase diagrams of such models at finite- 
N. An interesting gapless spin liquid discovered some time ago in the large-iV limit 
of the self-conjugate SU(N) Heisenberg model with biquadratic interactions to thwart 
dimerization is the 7r-ffux state [HU [17]. At N = oo, where projection to exactly N/2 
fermions per site does nothing to the mean field state, the 7r-flux state corresponds to 
the ground state of fermions (at half-filling) hopping on the square lattice while sensing 
a (spontaneously generated) fictitious magnetic flux of tt per elementary plaquette[16j. 
It supports gapless linearly dispersing Dirac fermion excitations about two nodes in the 
reduced Brillouin zone. Spin-1 excitations, which are bilinears of the Dirac fermions, are 
therefore gapless at the wavevectors spanning the Dirac nodes: (tt, tc), (jr, 0), (0, tt), and 
(0,0). For the case of the ordinary nearest-neighbor SU(2) Heisenberg antiferromagnet 
on the square lattice, Gutzwiller projecting this mean field wavefunction leads to a 
variational ground state with power law decay of staggered spin correlations (as there is 
no long range magnetic order) and an energy (l/2)(Sg(i)S^(i + 8)) ~ —0.319 per bond. 
The true ground state of the Heisenberg model is known to be Neel ordered, with a 
spin moment of about 60% of the classical value, and (1/2) (Sg(i)S^(i + 5)) ~ —0.3346 
per bond. Since the 7r-flux state is close to the true ground state, both energetically 
and in its display of (quasi) long-range antiferromagnetic correlations, we focus on the 
part of the phase diagram of the SU(N) Heisenberg model (with possible additional 
biquadratic interactions) that is close to that of the 7r-flux phase. More precisely, we 
investigate instabilities of the 7r-flux state towards various translational- and SU(N)- 
symmetry breaking orders. This point of view was advocated in earlier studies of the 
N = 2 case [371 fl~8] . and in more recent work by Ghaemi and Senthil |32j . 

We begin with the SU(N) Heisenberg model in the absence of biquadratic 
interactions, and compare the resulting variational phase diagram with that from a 
recent quantum Monte Carlo (QMC) study of the same model [34]. We then turn to 
the nature of the spin and dimer correlations of the £77(4) projected 7r-flux state and 
compare the correlation functions to results from recent analytical large- N studies [29] 
and QMC calculations [34]. Finally, we examine the variational phase diagram of the 
SU (4) and 577 (6) models with a biquadratic interaction added to thwart instabilities. 
For N = 6, we find a (small) window of parameters where the projected 7r-flux state 
appears to be stable towards Neel, spin Peierls and broken-C ordering. Our work 
thus hints at the existence of a stable SU(6) gapless spin liquid phase in a simple 
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(a) Dimer (Columnar) (b) Box (Plaquette) 

Figure 4. The two candidate spin Peierls ground states of the 2D SU(N) spin models 
explored in this paper. The thick lines indicate bonds with a larger singlet expectation 
value |<^(i)^(j)>|. 

two dimensional microscopic spin model. 

3.1. Phase Diagram of the SU(N) Heisenberg model 

The Hamiltonian of the 2D antiferromagnetic SU(N) Heisenberg model is 

H 2D = J2S^)S?(j), (15) 

where denotes nearest neighbor sites on the square lattice. As discussed above, the 
projected 7r-flux wavefunction with N/2 fermions at each site is an attractive starting 
point to describe the ground state of the model. The mean field ansatz for the 7r-flux 
phase is 

influx = E *ii {^ftfj* + h -c-) (16) 

(ij) 

with tij = t and a gauge choice of a iti+ x = |(— l) Xi +^ and a iii+ $ = — 1(— V) XiJrVi . Here 
we examine the instability of the variational state obtained by Gutzwiller projecting 
the ground state of this mean field Hamiltonian towards Neel and spin Peierls ordering. 
(We have also checked that there are no instabilities to time-reversal symmetry broken 
states or states with broken charge-conjugation symmetry. These wavefunctions have 
higher energy, so the ground states exhibit only Neel or spin Peierls order.) 

To account for the possibility of Neel ordering we modify the mean-field Hamiltonian 
with the addition of a SU(N) symmetry breaking perturbation that favors two-sublattice 
ordering, with any chosen set of N/2 flavors favored on one sublattice, and the remaining 
N/2 flavors on the other sublattice: 

H neel = H^-h N £ {-If^flfia+hH £ (-l) Xi+Vi flU -(17) 

i,a<N/2 i,a>N/2 
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Figure 5. Energy minimization plots for even N = 2 to 10. The x-axis is the 
variational parameter appropriate to the broken symmetry being studied. The y-axis 
is the (dimensionless) energy {{H 2 d) for the model of equation [15]) . We conclude that 
the SU(N) Heisenberg model exhibits Neel order for N = 2 and 4, and spin Peierls 
ordering of the columnar dimer type for N > 4. 



In order to study spin Peierls ordering, we focus on two different types of broken 
symmetry states, "dimer order" (more precisely, columnar dimer order) and "box order" 
(also called plaquette order). In the dimer state, spins prefer to form singlets on nearest 
neighbor bonds, and the bonds organize as shown in figure [U^a). Three other equivalent, 
but distinct, states are obtained by x— translations and n/2 rotations of the displayed 
pattern. The mean field ansatz for the preprojected wavefunction of the dimer state 
is obtained by modulating such that t^+z = 1 + <5d(— l) Xi and U s i+$ = 1. The box 
state has a different broken symmetry; the strength of singlet bonds is shown in figure 
Hfb). Three other equivalent box states are obtained by x— and y— translations of the 
displayed pattern. For the box state, we modulate Uj such that t^+i = 1 + <5b(— l) Xi 
and t iti+§ = 1 + 5 B {-l) y \ 
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Figure 6. Spin-spin correlation function of the 2D projected 7r-flux wavefunction for 
N = A and a system with L 2 sites. The label C(max) denotes the correlation function 
evaluated for points with the maximal separation, L/y/2, on the L x L square lattice. 



The mean field box and dimer states have a single particle gap to fermionic 
excitations and thus also a spin gap. These spin Peierls ordered states as well as the 
Neel state are invariant under a particle-hole transformation (followed by a global 577(4) 
spin rotation in the case of the Neel state), and thus are at half-filling. We project these 
mean field ansatz to obtain variational spin wavefunctions for the Heisenberg model. 
The phases of the SU (N) Heisenberg model are then obtained by looking for the state 
with the lowest variational energy. As summarized in figure we find that the Neel 
ordered state has the lowest energy for N = 2 and 4, while the dimer state (i.e., columnar 
dimer) state has the lowest energy for N > 4. For N = 2, this result is in agreement 
with other numerical work [52j [531 IM]- The presence of spin Peierls order for large 
values of N is in agreement with 1/N calculations [16]. The same pattern of (columnar) 
dimer order was also predicted for various representations of SU (N) antiferromagnets 
in large-N calculations by Read and Sachdev[5~4"l 155] ; these predictions have received 
some numerical support [56], with Neel order reported for N < 4 and dimer order for 
N > 4, identical to the phase diagram in figure [5j 

3.2. Correlation Functions of the Projected 577(4) tt-FIux State 

To make contact with a recent QMC study of the 577(4) Heisenberg model [54] , we turn 
now to the spin-spin and dimer-dimer correlations of the projected 7r-flux wavefunction. 
The analysis of these correlations is done in a manner similar to that in ID, except that 
we focus only on the strong nonzero-wavevector component (near (tt, it) for the spin 
order and near (tt, 0) for the dimer correlations) and ignore the uniform components. 
The uniform component of the spin-spin correlation in 2D is expected to decay quickly, 
as ~ 1/t 4 , and is therefore numerically harder to evaluate. 

The results for the finite size scaling of the staggered spin-spin correlation function 
are shown in figure El together with the correlation function results for the 12 x 12 
system. We find that the staggered spin-spin correlation function decays as l/r a with 
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a s (4) = 3.0(4). This value is in excellent agreement with large- iV calculations [2"9~| |2~T] 
that find a s (N) = 4—128/ (37r 2 iV). However, the exponent a s (4) is much larger than that 
found in a QMC simulation by Assaad Although the QMC calculation suggested 
a spin liquid ground state of the 7r-flux type for the SU(4) Heisenberg model, Assaad 
found a s ~ 1.12. The origin of the large difference needs further exploration, and we 
speculate on a possible reason for the discrepancy in the final section. 

A numerical analysis of the dimer-dimer correlation function at Q = (tt, 0) for the 
N = 4 projected 7r-flux wavefunction yields 0^(4) = 2.1(8); the larger error on the dimer- 
dimer correlation function exponent stems from having fewer Monte Carlo samplings of 
the wavefunction since the dimer-dimer correlation function takes more time to evaluate. 
We conclude that Gutzwiller projection strongly enhances both the spin-spin and the 
dimer-dimer correlations relative to the mean field result. In this sense, the projected ir- 
flux state is indeed the "mother of many competing orders" [29J! However, we have not 
confirmed yet that the projected 7r-flux wavefunction is an algebraic spin liquid phase 
with enlarged SU(2N) symmetry leading to a<f = a s [29J. We are currently carrying out 
further numerical calculations to reduce the error bars on the exponents and to better 
test this prediction quantitatively. 

3.3. Adding Biquadratic Interactions for N = 4 and N = 6 

The inclusion of the biquadratic interaction lead to a rich phase diagram in the case of 
the 577(4) spin chain. Motivated by this physics, we pursue here a variational study 
of the 2D square lattice model with Heisenberg bilinear and biquadratic interactions 
for SU(A) and SU(Q). As in ID, the 2D model is defined by the nearest-neighbor 
Hamiltonian: 

H 4 = cos£ £ S$(i)S£ti) + ^ £ WW^C*)] 2 , (18) 
(hi) (hi) 

where refer to nearest neighbor sites on the 2D square lattice. Exact 

diagonalization and the density matrix renormalization group are inadequate tools to 
study the phase diagram of the two dimensional model. For spin Hamiltonians without 
a sign problem, QMC has proved to be the most reliable numerical tool in 2D, but the 
case of a positive biquadratic spin interaction cannot be studied reliably because it is 
a frustrating interaction that introduces the sign problem. Given the success of the 
variational approach in ID we have reason to hope that the method may also provide 
a good guide to two dimensions, although we explore only a limited class of variational 
states. We consider the fully polarized ferromagnet TM., the (columnar) dimer state 
T>, the Neel antiferromagnet A/", the broken charge-conjugation symmetry state C, and 
the projected 7r-flux state n. 

T M. 1 Ferromagnet: The ground state of a fully polarized ferromagnet in 2D breaks 
global SU(N) spin symmetry but no lattice symmetries. An exact ground state 
wavefunction may be constructed by placing any two of the four fermion flavors on 
each lattice site, each site having the same two flavors. All other ground states can be 
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obtained by global SU(N) rotations of this state. Exactly as in ID, because the Pauli 
exclusion principle prevents any fermion hopping in the T ' M state, it is straightforward 
to show that (S|(i)S"£(i+l)) = N/A and ([S£(i)S£ (i + l)] 2 ) = N 2 /1Q, so that the exact 
ground state energy per site is 

N N 2 

e^(N) = -cos6 + — sm6. (19) 

J), Dimer: The columnar dimer state appeared in the phase diagram of the pure bilinear 
Heisenberg model and hence was discussed above in subsection 13.11 The only difference 
here is that we minimize the energy with respect to 5d at each value of 9. The energy 
of the perfectly dimerized state for N = 4, 6 is 

e™(5 D = l,iV = 4) = -^cos# + ^sin# (20) 

21 1197 

e™{5 D = 1, N = 6) = -— cos# + — — sinO. (21) 

4 oO 

A/", Neel: The Neel state with long range antiferromagnetic order was also discussed 
above. Here we minimize the energy with respect to the staggered magnetic field h^ at 
each 6. We emphasize that the classical antiferromagnetic state (obtained in the limit 
hjv — > 00) has a much higher energy than the optimal antiferromagnetic ground state 
at finite in the relevant region of the phase diagram. 

C, Charge-conjugation symmetry broken: The C phase is obtained by projecting the 
ground state of the mean field Hamiltonian 

H° c = H^ aux -t* J2 (-ir +m [fth + 2 ia + fl a fi + 2ya + h.C.}, (22) 
i,a=l..A 

the 2D generalization of equation [5j The intra-sublattice hopping t* gaps out the Dirac 
nodes of fermionic excitations in the 7r-flux state and thus leads to a spin gap in the 
mean field spectrum. 

n, The 7r-flux state: This gapless spin liquid state was also introduced earlier. This state 
does not have any variational parameters and is thus the most constrained state of the 
wavefunctions that we study. We leave for future work possible variational modifications 
of the state that preserve all lattice and spin symmetries. Since the 7r-flux state is a 
gapless state, it is important to study it under conditions such as particular system sizes 
that permit gapless nodes to appear at the mean field level. We find that the phase is 
artificially stabilized on lattices that do not permit the nodal wavevectors. 



3.4. Phase diagram of the SU(4) spin model 

We obtain the variational phase diagram shown in the left panel of figure [7] from 
an evaluation of the energy of the various 577(4) states. The ground state appears 
to generically exhibit broken symmetry. The T> — FAi, T> — M and the C — TM. 
transitions appear strongly first order due to level crossings. As in ID, the ground 
state is nearly completely dimerized at the T> — TM. transition. The location of the 
T> — TM. transition is therefore simply and reliably estimated by studying a dimer 
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Figure 7. Phase diagram of the SU(4) and SU(6) spin models of equation IT51 with 
Heisenberg bilinear and biquadratic interactions, on the 2D square lattice. Phases 
appearing here are the Neel phase, the ferromagnet, the columnar dimer phase, the 
broken charge-conjugation symmetry phase (C), and the 7r-flux spin- liquid phase (II). 
The thick line labeled "II?" in the N — 4 phase diagram indicates that the projected 
7r-flux state could either be stable over a thin sliver region (0.18 < 6/tt < 0.20) or 
instead there may be a direct transition at 8/ir sw 0.19 between the Neel and C phases. 
See text for details. 

product wavefunction as the variational state for T>. Setting ef^f = e™{o~D = 1) leads 
to ®d-tm ~ taa _1 (18/13) « -0.7tt. 

Due to limits on numerical accuracy and on system sizes (up to 10 x 10 for variational 
optimization), we are unable to determine whether there is a direct M — C transition 
at 6/n m 0.19 or a thin sliver of the 7r-nux phase that intervenes between these two 
phases for 0.18 < 6/tt < 0.20. Since the Neel and C states are both deformations of the 
7r-fiux state, it is possible for a direct continuous transition to occur between them, with 
the projected 7r-flux state being a possible candidate for the critical point. However 
the addition of variational parameter(s) to improve short distance correlations of the 
projected 7r-flux state could in fact stabilize this state. We are examining this issue 
more carefully, and comment further on this point in the final section. 

3.5. Phase diagram of the SU(6) spin model 

The variational phase diagram is shown in the right panel of figure [7J Strikingly the 
Neel phase is completely replaced by the dimerized phase at N = 6, and the biquadratic 
interaction appears to stabilize the 7r-flux state over a small window of 9. Of course 
within the variational approach one cannot rule out the possibility that II may be 
unstable to some other more complicated or exotic broken symmetries that have not 
been considered. 

The T> — FAi, T> — M and the C — TM. transitions appear strongly first order 
again due to level crossings. Since the ground state is nearly completely dimerized at 
the T> — TM. transition the location of this transition is reliably estimated by studying a 
dimer product wavefunction as the variational state for D. Setting e^ffi = e™(5[) = 1) 
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leads to B^_ FM « tarr^O/llO?) « -0.83vr. 
4. Summary and Discussion 

We have used Gutzwiller projected variational wavefunctions to deduce phase diagrams 
of SU(N) antiferromagnets with Heisenberg bilinear and biquadratic interactions in 
one and two spatial dimensions. In one dimension, the SU(4) variational phase diagram 
is in very good agreement with exact results. The spin and dimer correlations of the 
projected Fermi gas wavefunction with N fermion flavors are also in reasonably good 
agreement with 1/N calculations and exact results. Based on these results, the projected 
free fermion state with N fermion flavors appears to provides a good approximation of 
the critical points of SU(N) spin chains, and in particular it is a good description of 
the critical point between dimerized and broken charge-conjugation symmetry phases 
in the SU(A) model. 

On the two dimensional square lattice the pure bilinear Heisenberg model exhibits 
Neel order for N = 2 and 4 and columnar dimer order for N > A. Biquadratic 
interactions of positive sign appear to destabilize the Neel state, as the Neel order 
diminishes and gives way to a broken charge-conjugation symmetry phase via either a 
small sliver of the 7r-flux spin liquid or by a continuous transition that is well described 
by the projected 7r-flux state. 

The spin and dimer correlations of the projected 577(4) 7r-flux state are in 
reasonable agreement with analytical 1/N calculations. However the spin-spin 
correlations are quite different from those reported in QMC calculations by Assaad 
[34] . While that study finds a spin liquid ground state for the SU (4) Heisenberg model, 
apparently of the 7r-flux type, the spin correlations decay much more slowly than those 
predicted on the basis of 1/N calculations or our variational calculation. Based on the 
variational study of model equation [TS] with biquadratic interactions, we find that there 
could either be a thin sliver of the flux phase or a direct continuous Neel-C transition 
at#/7r~0.19. In the exact phase diagram this transition point, or the sliver of the flux 
phase, might occur even closer to the pure bilinear Heisenberg point. If this in fact is 
the continuous Neel-C transition with a 7r-flux state at the transition point (or 

a direct Neel-n transition if a region of stable 7r-flux spin liquid exists) could strongly 
influence the ground state of the pure S77(4) Heisenberg model as studied by QMC 
|34j . This hypothesis suggests that it may be numerically difficult to tell whether the 
correct ground state of the S77(4) Heisenberg model is a 7r-flux spin liquid or a Neel 
ground state with a much reduced staggered magnetization. It could also account for 
the discrepancy in the spin correlations between the QMC on one hand and the 1/N and 
variational calculations on the other. Further studies of the £77(4) Heisenberg model 
with biquadratic interactions might shed light on this issue. 

The 577(6) model does not appear to support a Neel phase at all. Instead 
biquadratic interactions open up a small window of 7r-flux phase between the dimerized 
and broken charge-conjugation symmetry phases. We checked for instabilities of 
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the 7r-flux state towards Neel order (characterized by non-zero (Sg(i))), dimer order 
(modulations in (TrS(i)S(j))) and C-breaking (modulations in (TrS(i)S(j)S(k))) and 
found it to be stable against all three. However, we cannot rule out instabilities 
towards other more exotic broken symmetries characterized by more complicated order 
parameters. While the dimerized ground state at Heisenberg point also appears to be 
close to a spin liquid phase from our phase diagram, it is less likely to be influenced by 
proximity to such a critical point as the dimerized phase has a spin gap rendering it 
more stable to critical fluctuations than the Neel phase. This picture is consistent with 
Assaad's QMC results, with a dimerized ground state reported for the SU(6) Heisenberg 
model. 

Finally, an exact C-breaking ground state of a 2D SU (8) spin model is known at a 
special point in parameter space [10] and it could be used as an additional test of the 
variational approach, which we have shown to be quite successful in describing a wide 
class of SU(N) antiferromagnets in one and two dimensions. 
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